require(ggplot2)
require(ggthemes)
require(tcltk)


rm(list=ls())


####  This script takes in the processed data from the force transducer and uses it to estimate stiffness of the fish.
    # First it log transforms the force data
    # Then it takes the log of the peak force necessary to bend the fish to each point and regresses that against the non log-transformed body curvature.
    # These slopes and R^2 values are recorded. They are corrected for size in the stats script, then used in the final analysis.


dir <- "Force data"
setwd(dir)


#Load in curvature data during 3-point bending tests
bend3PData <- read.csv("final3PBendCurvatureData.csv")


#Read the names of the .csv files in the current folder
fileNames <- list.files(pattern="*.csv", all.files=TRUE, no..=TRUE)
#Open all the files
list2env(
  lapply(setNames(fileNames, make.names(gsub("*.csv$", "", fileNames))), 
         read.csv,header=FALSE), envir = .GlobalEnv)


#Create data frame to hold data
finalDataStiffStickle <- data.frame(Individual=character(), stringsAsFactors = FALSE,
                                    sLslope <- double(),sLR2 <- double(),
                                    F3.5mm <- double(),Treatment=character(), Animal=double())



for(i in 1:length(fileNames)){
  nameNow <- fileNames[i]
  #Get rid of ".csv" in the filename
  nameNow2 <- substr(nameNow,1,nchar(nameNow)-4)
  #Add an x before the filename (because R data frame names can't start with a number)
  nameNow3 <- paste('X',nameNow2,sep='')
  #Open the file
  dataNow <- eval(parse(text= nameNow3))
  
  names(dataNow) <- c('Force','Time')
  
  
  #Extract video data from filename
  splitStrings <- strsplit(nameNow2,fixed=FALSE, split = '_')
  treatNow <- splitStrings[[1]][3]
  aNow <- as.integer(substr(splitStrings[[1]][4],2,3))
  
  
  #Create a variable that contains the displacement of the transducer
  dataNow$Displacement <- c(seq(0,3.5,.5),0)
  
  #correct the data by the force at 0 displacement.
  zeroSet <- dataNow[1,1]
  corrData <- dataNow[c(seq(2,8,1)),]
  corrData$Force <- dataNow[c(seq(2,8,1)),1]-zeroSet
  

  
  
  #find matching curvature data from 3-point tests
  p3CurveDF <- as.vector(bend3PData[which(bend3PData$Treatment==treatNow & bend3PData$Animal== aNow),])
  curvesNow <- as.vector(p3CurveDF[,c(seq(5,11,1))])
  
  
  corrData$Curvature <- t(curvesNow)
  
  
  #plot the graph of force vs. curvature
  ggplot(corrData, aes(x=Curvature ,y=Force)) + 
    geom_point(aes(),size=4,alpha=1,color='blue') +
    labs(x="Curvature (1/cm)", y="F (g)", title=c(nameNow2," Stiffness")) +
    theme_few()

  
  
  #Take log of force, then find linear slope of this corrected plot. (aka exponent of un-log data)
  #plot(logForce~corrData$Curvature, main=c(nameNow2,"semi-log"))
  logForce <- log(corrData$Force)
  lmTemp2 <- lm(logForce~corrData$Curvature)
  abline(lmTemp2)
  semiLogR2 <- summary(lmTemp2)$r.squared
  tempLmDf2 <- summary(lmTemp2)$coefficients
  semiLogSlope <- tempLmDf2[2,1]
  
  #plot the semi-log graph
  ggplot(corrData, aes(x=Curvature ,y=logForce)) + 
    geom_point(aes(),size=4,alpha=1,color='blue') +
    stat_smooth(method = "lm", col = "red") +
    labs(x="Curvature (1/cm)", y="Log(F)", title=c(nameNow2," Stiffness")) +
    theme_few()
  
  F3.5mm <- corrData$Force[7]
  
  finalDataStiffStickle[i,1] <- as.character(nameNow2)
  finalDataStiffStickle[i,2] <- semiLogSlope
  finalDataStiffStickle[i,3] <- semiLogR2
  finalDataStiffStickle[i,4] <- F3.5mm
  finalDataStiffStickle[i,5] <- treatNow
  finalDataStiffStickle[i,6] <- aNow
}

#Note, these slopes aren't size corrected! I will add the size correction in the final analysis script

names(finalDataStiffStickle) <- c('Individual','SLslope','SLR2','F3.5', 'Treatment', 'Animal')


fileSaveTo <- tclvalue(tkgetOpenFile(default=dir))
write.csv(finalDataStiffStickle, file=fileSaveTo)

